function eq = solveOptimalLSAP(etag,taug,taum,z0)

% Inputs:
% Solves for the optimal LSAP given z0
% Assume parameters are such that there is a unique equilibrium when lsap=0

% Returns the equilibrium with the optimal lsap

% Approach: Rather pedestrian
% Rewrite a bit the code of the solveEquilibrium to solve a class of problems as opposed to a single problem
%(assuming uniqueness)

% Underlying parameters(loaded through globals)
global g rho sigma PStar taub tauh k l tau

% This determines the accuracy: smaller is more accurate
lAcc = 0.0001;

assert(z0<=1);  % We are only interested in cases with productivity shocks are below baseline
eq.z0 =z0;
eq.etag=etag;
eq.taug=taug;
eq.taum =taum;
% Calculate the baseline equilibrium with lsap=0 and check condition for uniqueness
ph0 = exp(rho + g - sigma^2/tauh);
% Store all parameters for alter access
eq.params.g = g; eq.params.rho = rho; eq.params.sigma = sigma;
eq.params.tauh =tauh; eq.params.taub = taub;
eq.params.tau = tau;

assert(ph0<1); 
% Otherwise, households are sufficiently risk toelrant and lsaps are not needed
eq.ph0 = ph0;
zBar0 = fsolve(@(z) tau - k*l/z*(taub-tauh) -  sigma^2/(rho+g),1);
eq.zBar0 = zBar0;

zTilde0 = l/ph0;
eq.zTilde0 = zTilde0;

% Check condition for uniqueness
assert(zTilde0 < zBar0);

% My conjecture is that if this holds for lsap=0 it will also hold for
% higher lsaps. But I haven't proven this conjecture---if we see a
% counterexample the code will return an error

%%%%%%%%%% This is where the code departs from solveEquilibrium
pArr = linspace(ph0,1,100); % Normalized price array
% Unlike before, we will reuse this array for different lsaps, since lsaps
% increase the lower bound, we shouldn't lose any solutions by doing this
% (But reusing the arrays helps to speed up calculations)

sharpeActualArr = (g + rho - log(pArr))/sigma; % This is actually common across lsaps

% Calculate banks' implied balance sheets, and the implied risk tolerance, for each p
alphaArr = max(0,k*(1-l./(z0*pArr)));
tauBarArr = alphaArr*taub+(1-alphaArr)*tauh;

% Calculate the required Sharpe ratio. This is the main thing that changes
% for different lsaps. Calculate a baseline level and then multiply it by (1-lsap)

sharpeRequiredArr0 = sigma./tauBarArr; % For zero lsap

% Now calculate for each lsap in an array
% Keep going as long as government's utility keeps increasing
lsapArr = [0:lAcc:1];

Ug_lsapArr = 0*lsapArr - Inf;
p_lsapArr = 0*lsapArr - Inf;
ind_lsapArr = 0*lsapArr - Inf;
for i=1:length(lsapArr)
    lsap = lsapArr(i);
    % Check the boundary conditions without calling for fsolve (to speed
    % things up)
    ph = exp(rho+g - (1-lsap)*sigma^2/tauh);
    if z0*ph<l
        % Corner solution where bank is bankrupt
       p = ph;
    elseif ph>1 | tau - k*l/z0*(taub-tauh) >  (1-lsap)*sigma^2/(rho+g)
       p=1; 
       % Corner solution without demand recession
    else
   % Find the unique intersection
   [sol,numCross] = findCross((1-lsap)*sharpeRequiredArr0,sharpeActualArr,pArr);
   assert(numCross==1); % Intersection is supposed to be unique
   p = sol{1}.x;
%   ind = sol{1}.xind;
   end
   Ug = (1-exp(-rho))*log(p) - 1/2*exp(-rho)*sigma^2*(etag/taug*(1+lsap/etag)^2 + 1/taum*(1-lsap)^2); % Utility with 0 lsap
   % Set the arrays
   Ug_lsapArr(i) = Ug;
   p_lsapArr(i) = p;
   % Stop if we haven't improved
   if i>1 
       if Ug < Ug_lsapArr(i-1)
       lsapMax = lsapArr(i-1);
       pMax = p_lsapArr(i-1);
       UgMax = Ug_lsapArr(i-1);
       break; 
       % Break out of the for loop
       % We could calculate for all of the array, but this shouldn't be necessary
       end
   end 
end

eq.lsapMax = lsapMax;
eq.pMax = pMax;
eq.UgMax = UgMax;
% Return also the arrays for inspection
eq.lsapArr = lsapArr;
eq.p_lsapArr = p_lsapArr;
eq.Ug_lsapArr = Ug_lsapArr;

% Return also the baselinea arrays
eq.pArr = pArr;
eq.sharpeActualArr = sharpeActualArr;
eq.sharpeRequiredArr0  = sharpeRequiredArr0;

% Calculate other equilibrium variables corresponding to the max
alphaMax = max(0,k*(1-l/(pMax*z0)));
tauBarMax = alphaMax*taub+(1-alphaMax)*tauh;
rfMax= max(0,rho + g - (1-lsapMax)*sigma^2/tauBarMax); % Accommodates corner solutions
eq.alphaMax = alphaMax;
eq.tauBarMax = tauBarMax;
eq.rfMax = rfMax;

function [sol,numCross] = findCross(y1Arr,y2Arr,xArr)
    % To accommodate corner solutions, choose y1Arr and y2Arr so that
    % if y1Arr<y2Arr throughout then xArr(end) is the solution (and vice versa)
    
    % Find all crossings of two relations
    % y1(x1) and y2(x2)
    % defined on the same array, xArr
    % Requirement: The relations are continuous so every crossing will be interpreted as an equilibrium
    % Requirement: y1Arr(xArr)
    % Output:
    % sol{i}.y and sol{i}.x have the solutions corresponding to crossing i
    % sol{i}.xind has the x-index corresponding to the solution with crossing i
    % numCross: Number of equilibria      
    %
    % Corner cases: If there is no crossing, assume a corner solution (depending on y1<y2 or y1>y2) 
    numCross =0; 
    ind1Arr = []; 
    lastDif = 0; % Keeps track of the last difference, y1-y2, at the last step
    % newDif will keep track of the difference at this iteration
for ind=1:length(xArr)
    x = xArr(ind);
    y1 = y1Arr(ind);
    y2 = interp1(xArr,y2Arr,x,'linear');
    newDif = y1-y2;
    if isnan(y2)
        % No intersection
        continue;
    elseif newDif==0
        display('Crossing found!');
        numCross = numCross+1;
        sol{numCross}.y = y1;
        sol{numCross}.x = x;
        sol{numCross}.xind = ind;
    elseif newDif<0 & lastDif>0
        % Crossing. Implies an equilibrium
        display('Crossing found!');
        % Find the exact solution
        numCross = numCross+1;
        % Do a linear interpolation based on the magnitude of the differences
        newD = abs(newDif); lastD = abs(lastDif);
        sol{numCross}.y = y1Arr(ind-1)*newD/(lastD+newD)+y1Arr(ind)*lastD/(lastD+newD);
        sol{numCross}.x = xArr(ind-1)*newD/(lastD+newD)+xArr(ind)*lastD/(lastD+newD);
        sol{numCross}.xind = ind;
    elseif newDif>0 & lastDif<0
        % Crossing. Implies an equilibrium
        display('Crossing found!');
        numCross = numCross+1;
        newD = abs(newDif); lastD = abs(lastDif);
        sol{numCross}.y = y1Arr(ind-1)*newD/(lastD+newD)+y1Arr(ind)*lastD/(lastD+newD);
        sol{numCross}.x = xArr(ind-1)*newD/(lastD+newD)+xArr(ind)*lastD/(lastD+newD);
        sol{numCross}.xind = ind;
    end
    lastDif =newDif;
end  

    % If at the end of the loop numEq is zero, make some adjustment for corner solution
    if numCross ==0
        if lastDif<0
            % This means that y1<y2 throughout
            % Set the highest level of xArr as a corner solution 
            numCross=1;
            sol{numCross}.y = y1Arr(end);
            sol{numCross}.x = xArr(end);
            sol{numCross}.xind = length(xArr);
        elseif lastDif >0
            % This means that y1>y2 throughout
            % Set the lowest level of xArr as a corner solution 
            numCross= 1;
            sol{numCross}.y = y1Arr(1);
            sol{numCross}.x = xArr(1);
            sol{numCross}.xind = 1;
        end
    end

end

end